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Stochastic simulations are one of the cornerstones of the analysis of dynamical processes on 
complex networks, and are often the only accessible way to explore their behavior. The development 
of fast algorithms is paramount to allow large-scale simulations. The Gillespie algorithm can be used 
for fast simulation of stochastic processes, and variants of it have been applied to simulate dynamical 
processes on static networks. However, its adaptation to temporal networks remains non-trivial. 

We here present a temporal Gillespie algorithm that solves this problem. Our method is applicable 
to general Poisson (constant-rate) processes on temporal networks, stochastically exact, and up to 
multiple orders of magnitude faster than traditional simulation schemes based on rejection sampling. 

We also show how it can be extended to simulate non-Markovian processes. The algorithm is easily 
applicable in practice, and as an illustration we detail how to simulate both Poissonian and non- 
Mar kovian models of epidemic spreading. Namely, we provide pseudocode and its implementation 
in C++ for simulating the paradigmatic Susceptible-Infected-Susceptible and Susceptible-Infected- 
Recovered models and a Susceptible-Infected-Recovered model with non-constant recovery rates. 

For empirical networks, the temporal Gillespie algorithm is here typically from 10 to 100 times 
faster than rejection sampling. 


AUTHOR SUMMARY 


When studying how e.g. diseases spread in a pop¬ 
ulation, intermittent contacts taking place between 
individuals—through which the infection spreads—are 
best described by a time-varying network. This object 
captures both their complex structure and dynamics, 
which crucially affect spreading in the population. The 
dynamical process in question is then usually studied 
by simulating it on the time-varying network represent¬ 
ing the population. Such simulations are usually time- 
consuming, especially when they require exploration of 
different parameter values. We here show how to adapt 
an algorithm originally proposed in 1976 to simulate 
chemical reactions—the Gillespie algorithm—to speed up 
such simulations. Instead of checking at each time-step if 
each possible reaction takes place, as traditional rejection 
sampling algorithms do, the Gillespie algorithm deter¬ 
mines what reaction takes place next and at what time. 
This offers a substantial speed gain by doing away with 
the many rejected trials of the traditional methods, with 
the added benefit of giving stochastically exact results. 
In practice this new temporal Gillespie algorithm is tens 
to hundreds of times faster than the current state-of-the- 
art, opening up for thorough characterization of spread¬ 
ing phenomena and fast large-scale applications such as 
the simulation of city- or world-wide epidemics. 


* cvestergaard@gmail.com 


I. INTRODUCTION 

Networks have emerged as a natural description of 
complex systems and their dynamics [T], notably in the 
case of spreading phenomena, such as social contagion, 
rumor and information spreading, or epidemics 
The dynamics of contagion processes occurring on a net¬ 
work are usually complex, and analytical results are at¬ 
tainable only in special cases mia. Furthermore, such re¬ 
sults almost systematically involve approximations US]. 
Numerical studies based on stochastic simulations are 
therefore necessary, both to verify analytical approxima¬ 
tions, and to study the majority of cases for which no 
analytical results exist. The development of fast algo¬ 
rithms is thus important for the characterization of con¬ 
tagion phenomena, and for large-scale applications such 
as simulations of world-wide epidemics nig. 

The Doob-Gillespie algorithm [BHTT] (also known as 
Gillespie’s Stochastic Simulation Algorithm—SSA or 
Gillespiedirect method), originally proposed by David 
Kendall in 1950 for simulating birth-death processes and 
made popular by Daniel Gillespie in 1976 for the simula¬ 
tion of coupled chemical reactions, offers an elegant way 
to speed up such simulations by doing away with the 
many rejected trials of traditional Monte Carlo methods. 
Instead of checking at each time-step if each possible re¬ 
action takes place, as rejection sampling algorithms do, 
the Gillespie algorithm draws directly the time elapsed 
until the next reaction takes place and what reaction 
takes place at that time. It is readily adapted to the 
simulation of Poisson processes on static networks p!2UT6] 
and can be generalized to non-Markovian processes HU. 

Systems in which spreading processes take place, e.g., 
social, technological, infrastructural, or ecological sys¬ 
tems, are not static though. Individuals create and break 
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contacts at time-scales comparable to the time-scales of 
such processes [TsHlO] , and the dynamics of the networks 
themselves thus profoundly affect dynamical processes 
taking place on top of them |2lH27]. This means that 
one needs to take the network’s dynamics into account, 
e.g., by representing it as a time-varying network (also 
known as a time-varying graph, temporal network, or 
dynamical network) [28]. The dynamical nature of time- 
varying networks makes the adaptation of the Gillespie 
algorithm to such systems non-trivial. 

The main difficulty in adapting the Gillespie algorithm 
to time-varying networks is taking into account the vari¬ 
ation of the set of possible transitions and of their rates 
at each time step. We show that by normalizing time 
by the instantaneous cumulative transition rate, we can 
construct a temporal Gillespie algorithm that is applica¬ 
ble to Poisson (constant rate) processes on time-varying 
networks. We give pseudocode and G++ implementa¬ 
tions for its application to simulate the paradigmatic 
Susceptible-Infected-Susceptible (SIS) and Susceptible- 
Infected-Recovered (SIR) models of epidemic spreading, 
for both homogeneous and heterogeneous [29| popula¬ 
tions. We verify the accuracy of the temporal Gillespie 
algorithm numerically by comparison with a classical re¬ 
jection sampling algorithm, and we show that it is up to 
^ 500 times faster for the processes and the parameter 
ranges investigated here. 

While Poissonian models are of widespread use, real 
contagion phenomena show memory effects, i.e., they are 
non-Markovian. Notably, for realistic infectious diseases, 
the rate at which an infected individual recovers is not 
constant over time [301EB. Social contagion may also 
show memory effects, e.g., one may be more (or less) 
prone to adopt an idea the more times one has been ex¬ 
posed to it. To treat this larger class of models, we show 
how the temporal Gillespie algorithm can be extended 
to non-Markovian processes. We give in particular an 
algorithm for simulating SIR epidemic models with non¬ 
constant recovery rates. 


II. RESULTS 


The following subsections present the main results of 
the article. Section EA] defines the stochastic processes 
which can be simulated using the temporal Gillespie al¬ 
gorithm, and describes the class of compartmental mod¬ 
els for contagion phenomena on networks—the class we 
will use in examples throughout this paper. Section [II B| 
gives a quick overview of the traditional rejection sam¬ 
pling algorithms. Section [TrC] outlines a derivation of the 
static Gillespie algorithm. Section m derives the tem¬ 
poral Gillespie algorithm for Poisson (constant-rate) pro¬ 
cesses. In Section [ilEj we validate the temporal Gillespie 
algorithm through numerical comparison with a rejec¬ 
tion sampling algorithm; we also compare their speeds for 
simulating SIR and SIS processes on both synthetic and 
empirical time-varying networks. Section [II Fj shows how 


the temporal Gillespie algorithm can be extended to sim¬ 
ulate non-Mar kovian processes; the approach is verified 
numerically and the speed of the non-Markovian tempo¬ 
ral Gillespie algorithm is compared to rejection sampling. 

Tables listing the notation used in the manuscript, de¬ 
tails on how Monte Carlo simulations were performed, 
and pseudocode for application of the temporal Gillespie 
algorithm are given in the Methods section. 


A. Stochastic processes on time-varying networks 

We define in this section the type of stochastic pro¬ 
cesses for which the temporal Gillespie algorithm can be 
applied. At the time of writing, the main domain of ap¬ 
plication of the algorithm is the class of compartmental 
models for contagion processes on time-varying networks, 
which we introduce below. For definiteness, algorithms 
detailing the application of the temporal Gillespie algo¬ 
rithm will concern this class of stochastic processes. 

In general, we consider a system whose dynamics can 
be described by a set of stochastic transition events. We 
assume that the system can be modeled at any point in 
time by a set, f^(t), of M{t) independent stochastic pro¬ 
cesses m, which we term transition processes; the rate at 
which the transition m takes place is denoted A^. The 
set 0(t) thus defines the possible transition events at time 
t and in general changes over time, depending on both 
external factors and the evolution of the system itself; the 
number of possible transitions, M(t), thus also generally 
changes over time, while Xm may or may not vary over 
time. For the classic “static” Gillespie algorithm to be 
applicable, U(t) is allowed to change only when a tran¬ 
sition (or chemical reaction in the context of Gillespie’s 
original article) takes place. For processes taking place on 
time-varying networks, the medium of the process—the 
network—also changes with time. As these changes may 
allow or forbid transitions, Cl{t) is not only modified by 
every reaction, but also by every change in the network. 
This is the domain of the temporal Gillespie algorithm, 
which only requires that the number of points in which 
0(t) changes be finite over a finite time-interval [32]. 

The assumption that the transition processes are in¬ 
dependent is essential to the validity of the Gillespie al¬ 
gorithm, as it allows the calculation of the distribution 
of waiting times between consecutive transitions. This 
assumption is not overly restrictive, as it only requires 
a transition process to be independent of the evolution 
of the other simultaneous transition processes. A tran¬ 
sition process may depend on all earlier transitions, and 
the current and past states of all nodes. As such, Gille¬ 
spie algorithms can notably be applied to models of co¬ 
operative infections and other non-linear processes such 
as threshold models and has even been applied to 
model the dynamics of ant battles [33] . 

Compartmental models of eontagion. In a network- 
based description of the population in which a conta¬ 
gion process takes place, an individual is modeled as 
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a node i (Fig. [^). A contact between two individuals 
taking place at time t is represented by an edge (i, 
in the graph describing the population at the instant t 
(Fig. [^). In a compartmental model, each node is in 
a certain state, which belongs to a fixed, finite set of q 
different states (compartments) [3]. A random variable 
Xi{t) G {Ai, A 2 ,..., Ag} specifies the state of the node i 
at time t (i.e. to which compartment it belongs). Nodes 
may stochastically transition between states, governed 
by the set of transition processes. One is usually in¬ 
terested in the evolution of the number of nodes in each 
state, which we denote Xi, X 2 ,..., Xq. 

As an example, we consider the classic SIR model of 
epidemic spreading with constant transition rates in a 
homogeneous population (rates are the same for all in¬ 
dividuals) (Fig. m- Here nodes can be in one of three 
states: susceptible, infected, and recovered, {tS,X, 7^}. 
Two different transition types let nodes change state: (i) 
a node i in the S state switches to the X state with rate 
kx{t)P (an S ^ X reaction), where kx{t) is the num¬ 
ber of contacts i has with nodes in the X state at time t 
(Fig. [^); (ii) a node i in the X state switches to the 7Z 
state at rate /i (an X ^ IZ reaction). 

In general, the transition processes can be divided into 
three types: 


a. spontaneous transitions, which only depend on the 
current state of the node, Xi{t) (Fig. [ip)— e.g. an 
infected node recovers spontaneously in the SIR 
model (Fig. [^); 

b. contact-dependent transitions, which may take 
place only when the node i is in contact with other 
nodes in a given state; it thus depends on the 
states Xj of the nodes j currently in contact with i 
(Fig. [^)— e.g. a susceptible node may be infected 
upon contact with an infectious node in the SIR 
model (Fig. [^). 

c. mixed transitions, which take place spontaneously, 
but may depend on the node’s past and current 
contacts (Fig. [^)— e.g. in rumor spreading, an 
individual may learn on his own that a rumor is 
false (spontaneous) or may be convinced by another 
individual who knows the rumor is false (contact- 
dependent). 


This division is important for practical application of 
the temporal Gillespie algorithm as transition processes 
of type a need only be updated after a transition has 
taken place, and processes of type c need only be updated 
whenever a relevant contact takes place, but not at each 
time-step. Using this distinction is crucial to obtaining 
the large speed-increase that the temporal Gillespie algo¬ 
rithm offers over rejection sampling, as discussed below 
(Secs. 


IID 


and 


HE). 



FIG. 1. Schematic representation of a compartmental 
contagion process on a network. (A) Illustration of a 
contagion process evolving on a time-varying network. Nodes’ 
colors correspond to their current state; edges denote current 
contacts between nodes; edge colors correspond to: black: no 
contagion may take place over the edge, red: contagion takes 
place during the present time-step, and red-to-blue gradient: 
contagion is possible but does not take place. (B) Example: 
reaction types in the SIR model. (C) Spontaneous reaction: 
a node i may spontaneously transition from its current state 
Xi to x'i with rate Am- (D) Contact-dependent reaction: a 
node i may transition from its current state Xi to x'^ with 
rate Am upon contact with the node j in state Xj . (E) Mixed 
transition: a node i may spontaneously transition from its 
current state Xi to another state, x- with rate Am; contact 
with another node j, in state Xj, may alter the transition rate 
of m, Am A(^. After the contact ends, the transition 

rate may revert to Am, remain unchanged, or change to a 
third value. 
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FIG. 2. Example of how the temporal Gillespie algorithm works for an SIR process on a time-varying network. 

We consider the time-varying network of Fig[^ (Network)—time evolves along the vertical axis; a rejection sampling algorithm 
considers each transition process at each time-step individually (Transitions); the temporal Gillespie algorithm considers the 
integrated cumulative transition rate of all transition processes, L(t;0), and compares it with the random normalized waiting 
times r/ (Normalized time). A transition takes place whenever L(t; 0) = ^^ 1 = 0^1 foi* any g G N. The temporal Gillespie 
algorithm works as follows. (A) The first normalized waiting time r[ is drawn from an exponential distribution with unit rate 
[r[ ~ Exp(l)] (Normalized time). From the state of the network at the first time-step, the set of possible transitions Q(0) is 
found (Transitions), and from this the cumulative transition rate A(0) is calculated. The integrated cumulative transition rate, 
L(At; 0) = A(0)At is compared to r[. If, as in the present example, A(0)At < r[ the algorithm is advanced to the next time- 
step. (B) and (G) The set of possible transitions Q(tn) and the cumulative transition rate A{tn) is updated at each following 
time-step n; if L(tn; 0) = At A{ti) is still smaller than , the algorithm is advanced to the next time-step. (D) During the 

first time-step n** where L(tn**; 0) > , a transition takes place. The exact time of this transition, t**, is given by Eq. (12) and 

the transition m that takes place is chosen among the possible transitions in the given time-step with probability proportional 
to its transition rate Am [Transitions and Eq. ( |10| ]. (E) The transition changes the system (Network) and consequently the set 
of possible transitions, Q(t**), (Transitions); thus Q(t**) and A(t**) must be updated, which in turn changes the remainder of 
L(tn**;0) (Normalized time). A new normalized waiting time is then drawn, T 2 ~ Exp(l); if L(tn**;0) < + T 2 , no further 

transitions takes place during the time-step and the algorithm is advanced to the next time-step (note that multiple transitions 
may occur during the same time-step). (E) The above procedure is reiterated. 
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B. Rejection sampling for Monte Carlo simulations 

A straightforward way to simulate a stochastic process 
is to use a rejectiou sampliug algorithm, akiu to the clas¬ 
sical Metropolis algorithm. Here oue divides the time- 
axis iu small time-steps At, where At should be choseu 
sufficieutly small such that this discretizatiou does uot iu- 
flueuce the outcome of the process siguificautly; ou time- 
varyiug uetworks, the choice of At ofteu comes uaturally 
as the time-resolutiou at which the uetwork is measured 
or simulated (Fig. 

At each time-step t = 0, At, 2At,..., we test whether 
each possible trausitiou m G Q{t) takes place or uot. Iu 
practice this is doue by drawiug a raudom uumber 
that is uuiformly distributed ou [0,1) for each m aud 
compariug it to Xm^t: if < A^At the reactiou takes 
place, if Vm > Xm^Xt uothiug happeus [Fig. (Trausi- 
tious)]. (Note that oue should techuically compare to 
1 — exp(A^At) to eusure that A^ defiues a proper trausi¬ 
tiou rate for fiuite At. However, the two procedures are 
equivaleut iu the limit At ^ 0.) 

From the desigu of the rejectiou sampliug algorithm we 
see that the proportiou of trials that are rejected is equal 
to a weighted average over {1 — A^At}^. Thus, siuce we 
require A^At ^ 1 iu order to avoid discretizatiou errors, 
the vast majority of trials are rejected aud the rejectiou 
sampliug algorithm is computatioually iuefficieut. 


C. Gillespie algorithm on static networks 

The Gillespie algorithm lets us perform stochastically 
exact Moute Carlo simulatious without haviug to reject 
trials. For Poissou processes ou static uetworks, it works 
by recoguiziug that the waitiug time betweeu two cousec- 
utive trausitious is expoueutially distributed, aud that 
each trausitiou happeus with a probability that is pro- 
portioual to its rate. 

Specifically, the (survival) probability that the trausi¬ 
tiou m has uot takeu place after a time r siuce the last 
trausitiou eveut is 


5™(t) = . (1) 

Siuce each trausitiou takes place iudepeudeutly, the prob¬ 
ability that uo eveut takes place duriug the iuterval r 
siuce the last eveut is 


beiug the uext reactiou that takes place aud that it takes 
place after exactly time r is equal to Pm{^) = 

The static Gillespie algorithm thus cousists iu drawiug 
the waitiug time r ^ Exp (A) uutil the uext trausitiou 
aud theu drawiug which trausitiou m takes place with 
probability tt^ = A^/A. [Here r ^ Exp(A) is short for: 
r is expoueutially distributed with rate A.] 


D. Temporal Gillespie algorithm 


Eor processes takiug place ou time-varyiug uetworks 
however, the set of trausitiou process, chauges with 
time iudepeudeutly of the trausitiou eveuts, e.g., for the 
case of au SIR process uodes may become iufected ouly 
wheu iu coutact with au iufected iudividual (Eig. 

This meaus that the survival probability does uot reduce 
to a simple expoueutial as iu Eq. 0 ; it is iustead giveu 
by 


Sm{T;t*) = exp (- 


Im XrfJidt 


( 3 ) 


where t* is the time at which the last trausitiou took 
place, t** = t* + r is the time wheu the uext trausitiou 
takes place, aud lm{t) is au iudicator fuuctiou that is 
equal to oue wheu the process m may take place, e.g., 
wheu two giveu uodes are iu coutact, aud zero wheu m 
may uot take place. The meauiug of Im is exemplified 
iu Eig. [T}\: the uode i may be iufected by the iufectious 
uode j ouly wheu the two uodes are iu coutact; if we let 
m deuote this trausitiou process, lm{t) is theu oue for 
t = At, 3At, 4At aud zero for t = 0, 2At. 

Note that for processes takiug place ou adaptive time- 
varyiug uetworks, whose chauges ouly depeud ou the pro¬ 
cess itself, I^(t) ouly chauges wheu a trausitiou takes 
place aud Eq. Q reduces to Eq. 0. This meaus that 
from the poiut of view of the algorithm, such uetworks 
are effectively static aud the classic “static” Gillespie al¬ 
gorithm may simply be used there [nnn]. 

We uow cousider the geueral case where Q{t) may 
chauge iudepeudeutly of the proc esses evolviug ou the 
uetwork (as described iu Sec. HA). Usiug, as iu the pre¬ 
vious sectiou, that trausitiou processes are iudepeudeut, 
we cau write the probability that uo eveut takes place 
duriug au iuterval r (the waitiug time survival fuuctiou): 


- (2) 

m 

where A = Ylm=i is the cumulative trausitiou rate. 
The above result is obtaiued by usiug the fact that while 
Q aud M do depeud ou t, they ouly chauge wheu au eveut 
takes place aud uot iu-betweeu. They cau thus be treated 
as coustaut for the purpose of calculatiug the waitiug 
time betweeu eveuts. The distributiou of the waitiug 
times r is theu giveu by the probability deusity p{r) = 
Ae”^'^, while the probability deusity for the reactiou m 


S{T-,t*)= n SUr-,n 

= exp|-y; f I^{t)Xmdt\ , ( 4 ) 

\ mEQ ^ / 

where Q deuotes the set of all possible trausitious (trau¬ 
sitiou processes) ou the iuterval betweeu two trausi¬ 
tiou eveuts, (t*,t**], i.e., ft is the uuiou over ft{t) for 
t G (t*, t**], aud M is the total uumber of trausitiou pro¬ 
cesses ou the same iuterval (the size of ft). We switch 
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the sum and the integral in Eq. @ to obtain 

S{T-,t*) = exp (- f . (5) 

V men J 


the remainder of the n*th time-step left after the last 
transition; the second term is equal to L(tn+i; tn*+i). A 
new transition takes place during the time-step n** where 
> r' (Fig. [^); the precise time of this new 

transition is 


Finally, using that lm{t) = 0 for all m ^ we may 

write 


S(t; t*) = exp 



where 

A(i) = Am 

men(t) 


(6) 

(7) 


is the cumulative transition rate at time t. 

The dynamics of empirical time-varying networks is 
highly intermittent and we cannot describe analyti¬ 
cally. This means that we cannot perform the integral of 
Eq. ^ to find the waiting time distirbution directly. We 
may instead normalize time by the instantaneous cumula¬ 
tive transition rate, A(t): We define a unitless normalized 
waiting time between two consecutive transitions, r', as 


t' = L (t**;t*) = A{t)dt , (8) 

i.e., equal to the cumulative transition rate integrated 
over (t*, t**]. The survival function of r' has the following 
simple form: 


S'(t') = exp(-T') . (9) 

The time t** when a new transition takes place is given 
implicitly by L(t**; t*) = r', while the probability that m 
is the transition that takes place at time t = t** is given 
by: 

= lm{t)Xm/A{t) . (10) 


This lets us define a Gillespie-type algorithm for time- 
varying networks by first drawing a normalized waiting 
time r' until the next event from a standard exponen¬ 
tial distribution [i.e. with unit rate, r' ^ Exp(l)], and 
second, solving L(t;t*) = r' numerically to find t**. In 
practice, since A(t) only changes when a transition takes 
place or at tn = nAt with n G N, we need only compare 
r' to 


Utn+i;t*) = {tn^+l-t*)Ht*) + At Hti) AW 

i=n*-|-l 

for each time-step n (Fig. -C). Here n* is the time- 
step during which the last transition took place, and 
A(t*) is the cumulative transition rate at t*, immediately 
after the last transition has taken place. The first term of 
Eq. (0 is the cumulative transition rate integrated over 


tin,** 





( 12 ) 


the reaction m that takes place is drawn with probability 
given by Eq. (Fig. pD ). We then update Q and A 
to f](t**) and A{t**) (Fig7[^), draw a new waiting time, 
r' ^ Exp(l), and reiterate the above procedure (Fig.[^). 

The algorithm can be implemented for conta¬ 
gion processes on time-varying networks as follows 
(see Methods for pseudocode for specific conta¬ 
gion models and https: //github. com/CLVestergaard/ 
TemporalGillespieAlgorithm for implementation in 
C++): 


1. Draw a normalized waiting time until the first 
event from a standard exponential distribution, 
r' Exp(l) (Fig. [^). 

2. At each time-step tn = nAt, with n = 0,1, 2,..., 
let O = ^{tn) and A = A(t^); here, only 


contact-dependent pr ocess es (type b. Sec. II A) and 
mixed (type c. Sec. IIA) processes that depend 
on contacts taking place at tn or tn-i need to be 
updated—an important point, as it lets the tempo¬ 
ral Gillespie algorithm be much faste r tha n rejec¬ 
tion sampling (see discussion in Sec. HE). Then, 
compare r' to AAt: 


if AAt < r'l Subtract AAt from r', continue to 
next time-step and repeat 2 (Figs. 

IMl. 


if AAt > t': Let the reaction m take place, chosen 
from D with probability tt^ = A^/A. The 
fraction that is left of the time-step when the 
transition takes place is ^ = 1 — t'/{AA t) 
and the precise time of the transition is t** = 
tn +t'/A (Figs.j^ and[^). Next, update ft 
and A (Fig. [^); this time all transition pro¬ 
cesses should be up dated , as spontaneous pro¬ 
cesses (type a. Sec. II A) may change, emerge, 
or disappear when a transition takes place. 
Then: 


(a) draw a new normalized waiting time, r' ^ 

Exp(l) (Fig. If’); 

(b) compare r' to ^AAt: 

if r' > ^AAt: subtract ^AAt from r', 
continue to the next time-step and re¬ 
peat 2 (Fig. W)- 

if r' < ^AAt: Another transition takes 
place during the present time-step (at 
time t*** = t** + r'/A, where t** is 
the time of the last transition during 
the same time-step): choose m from 
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ft with probability tt^ = A^/A; let 
^ ^ ^ — r'/AAt, and update and 
A. Repeat a) and b). 

By construction, the above procedure produces re¬ 
alizations of a stochastic process for which the wait¬ 
ing times for each transition follow exactly their cor¬ 
rect distributions. The temporal Gillespie algorithm is 
thus what we term stochastically exact: all distributions 
and moments of a stochastic process evolving on a time- 
varying network obtained through Monte Carlo simu¬ 
lations converge to their exact values. Rejection based 
sampling algorithms are stochastically exact only in the 
limit A^At —> 0. 

A large literature exists on the related problem of 
simulating coupled chemical reactions under externally 
changing conditions (e.g., time-varying temperature or 
volume) [35H4Q] . Most of these methods consider only 
external perturbations that can be described by an an¬ 
alytical expression. In this case the problem reduces to 
that of defining a static, yet non-Markovian, algorithm. 
Some methods, and notably the modified next reaction 
method developed by Anderson can be adapted to a 
completely general form of the external driving and thus, 
in principle, to simulate dynamical processes taking place 
on time-varying networks. These methods are based on 
a scheme that is conceptually similar to Gillespie’s direct 
algorithm, the next reaction method, proposed by Gib¬ 
son and Bruck [35]. The next reaction method draws a 
waiting time for each reaction individually and chooses 
the next reaction that happens as the one with the short¬ 
est corresponding waiting time. It then updates the re¬ 
maining waiting times, draws new waiting times (if ap¬ 
plicable), and reiterates. To generalize the next reaction 
method to processes with non-exponential waiting times, 
Anderson introduced the concept of the internal time for 
each transition process m- In the notation used in the 
present article it is defined as Tm{t) = and 

is thus equivalent to the normalized time, L(t,0), only 
for an individual transition process. 

By construction, the next reaction method needs to 
draw only one random number per transition event, 
where the Gillespie algorithm draws two. However, this 
reduction in the number of required random variables 
comes at a price: one must draw a random number for 
each individual transition process and keep track of, com¬ 
pare, and update each of the individual waiting times. 
For chemical reactions, where the number of different 
chemical reactions is small (it scales with the number of 
chemical species), this tradeoff favors the next reaction 
method. However, for contagion processes on networks, 
each individual is unique (if not intrinsically, at least due 
to its position in the network). The number transition 
processes thus scales with the number of nodes and con¬ 
tacts, which favors the Gillespie algorithm as it does not 
need to keep track of each of them individually m- 

On time-varying networks (or for time-varying exter¬ 
nal driving) one must furthermore update relevant in¬ 


ternal times each time the network structure (external 
conditions) changes in the next reaction method. Ghem- 
ically reacting systems are usually close to being adia¬ 
batic, i.e., the external driving changes slowly compared 
to the time-scales of chemical reactions. Thus, the ad¬ 
ditional overhead related to updating individual internal 
times is practically negligible. However, the dynamics of 
temporal networks is highly intermittent and the time- 
scale of network change is typically smaller than the time- 
scales of relevant dynamical processes. Here one must 
thus update the internal times many times between each 
transition event, inducing a substantial overhead. Since 
the temporal Gillespie algorithm operates with a single 
global normalized waiting time, it handles these updates 
more efficiently. 

Finally, the modified next reaction method may in 
principle be extended to non-Markovian processe s tak ing 
place on time-varying networks (as treated in Sec. H F us¬ 


ing the temporal Gillespie algorithm). However, such an 
approach would, for each single transition, require solving 
numerically Eq. (13) of [37] for the internal waiting time 
of each individual transition process, taking into account 
the time-varying network structure, finding the shortest 
corresponding waiting time in real time, and then updat¬ 
ing the internal waiting times of all the other reactions, 
rendering the next reaction method even more inefficient 
in this general case. 


E. Comparison of Gillespie and rejection sampling 
algorithms 

Numerical validation. We compare the outcome of 
SIR and SIS processes on activity-driven time-varying 
networks m simulated using the temporal Gillespie al¬ 
gorithm to the outcome of simulations using traditional 
rejection sampling. For sufficiently small A^At, the out¬ 
comes are indistinguishable (Fig. see also Supplemen¬ 
tary Fig. 12 for an empirical network of face-to-face con¬ 
tacts in a high school), confirming the validity of the tem¬ 
poral Gillespie algorithm. Note that rejection sampling 
is only expected to be accurate for A^At ^ 1, while the 
temporal Gillespie algorithm is stochastically exact for 
all A^At; the results of the two algorithms thus differ 
when the assumption A^At ^ 1 does not hold (Supple¬ 
mentary Fig. [^. 

Comparison of simulation speed. Next, we compare 
the speeds of the temporal Gillespie and the rejection 
sampling algorithms for SIR and SIS processes (see Meth¬ 
ods for details on how simulations were performed). Fig¬ 
ure shows that the temporal Gillespie algorithm is up 
to multiple orders of magnitude faster than traditional 
rejection sampling. These results are confirmed by simu¬ 
lations on empirical time-varying networks of face-to-face 
contacts (Fig. Table |^. The speed gain is higher for 
larger systems Compare N = 1 000 to N = 100 in Fig.[^ 
We also see that the speed gain is larger the sparser the 
network is. This is because the calculation of the contacts 







FIG. 3. Comparison of numerical results from temporal Gillespie and rejection sampling algorithms. (A) Mean 
number of nodes in each state of the SIR model as function of time. (B) Distribution of final epidemic sizes (number of recovered 
nodes when / = 0) in the SIR model. (C) Mean number of nodes in each state of the SIS model as function of time. (D) 
Distribution of the number of infected nodes in the stationary state (t ^ oo) of the SIS model. All simulations were performed 
1 000 000 times with the root node chosen at random on an activity driven network m consisting of A" = 100 nodes, with 
activities ai = r]Zi, where 77 = 0.1 and Zi ~ for Zi G [0.03,1), and a node formed two contacts when active. Parameters of 

the epidemic processes were f3At = 10“^ and /xAt = 10“^. 


TABLE I. Summary statistics for empirical face-to-face 
contact networks from the SocioPatterns collabora¬ 
tion [45| : social setting; number of nodes in the network, 
N ; tota l dur ation of measurements, T ; average instantaneous 
degree, k{t). 


Setting 

AT 

T 

k(t) 

Reference 

Workplace 

92 

11 days 

0.004 

m 

Hospital 

80 

4 days 

0.064 

m 

High school 

327 

4 days 

0.063 

m 

Gonference 

399 

32 hours 

0.070 

m 


between susceptible and infected nodes at each time-step, 
necessary to determine the possible 5 —> X transitions, is 
the performance limiting step of the temporal Gillespie 
algorithm (see below). In the extreme case of a conta¬ 
gion model where all transitions are contact-dependent 
(type b, Sec. IIIAI ), such as the classic Maki-Thompson 
model of rumor spreading [42], the temporal Gillespie 
algorithm is approximately a factor two faster than the 
rejection sampling algorithm. 


Expected time complexity of the algorithms. We may 
gain insigth into the performance of the algorithms by 
considering their time-complexity, i.e., how their run¬ 
ning time scales with the input parameters of the sim¬ 
ulated system. Since the algorithms are used for Monte 
Carlo simulations, it is most interesting to consider the 
expected complexity given a set of parameters, i.e., the 
mean running time of an algorithm averaged over an 
ensemble of simulations, not the worst-case complexity 
which is usually considered for deterministic algorithms. 

The expected running time of the rejection sampling 
algorithm scales as 

grs = o (m O (M(t) 

^simu ) , (13) 

where O (x) denotes a term that is of order x, E(t) = 
Nk(t)/2 is the mean number of contacts per time-step, 
M(t) is the mean number of possible transitions at any 
instant, and rigimu is the number of time-steps simulated. 
For comparison, the expected running time of the tem¬ 
poral Gillespie algorithm is given by 

0TGA = O (m ^simu )+o(Q(t) 

^simu ) , ( 14 ) 
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FIG. 4. Comparison of the speed of the temporal Gillespie and the rejection sampling algorithms. Ratio between 
computational times 0rs and 0tga per single realization of a spreading process using rejection sampling and the temporal 
Gillespie algorithm, respectively: (A) for a SIR process and (B) for a SIS process on networks with different mean degree, 
k{t) [for empirical contact networks, k(t) ~ 0.004-0.07 (Table |^]. Networks consisted of A = 100 or A = 1000 nodes, with 
activities ai — r]Zi and Zi ~ for Zi G [0.03,1); a node formed two contacts each time it was active. For At = 20 s (as 

for the empirical data), /xAt 3 • 10“^ corresponds to a recovery time of roughly one week, typical of flu-like diseases. The 
infection rate was [3 = 10^/i for networks with k(t) = 0.002, [3 = 10^/i for networks with k(t) = 0.02, and [3 = 10/i for networks 
with k{t) = 0.2. (Details on how simulations were performed are found in Methods.) 


where Q{t) is the mean number of transitions that take 
place per time-step. 


The first term of the r.h.s. of Eqs. (13) and (14) 
correspond to the time needed for looking through the 
set of contacts at each time-step to determine the set 
of possible infections and are thus similar for the rejec¬ 
tion sampling and temporal Gillespie algorithms (with 
the temporal Gillespie algorithm incurring a small ad¬ 
ditional overhead related to calculating the cumulative 
transition rate and keeping track of of the normalized 
waiting time left till the next transition). For rejection 
sampling [Eq. (p!^], the second term corresponds to the 
determination of whether each of the possible transitions 
takes place at each time-step; for the temporal Gille¬ 
spie algorithm [Eq. (p^], the second term corresponds to 
drawing inter-event waiting times and which transitions 
that ta ke pla ce. Eor the S IR and SIS pr ocesses co nsidered 
above, M{t) = Ms^i(t) + /(t), where Ms^i(t) is mean 
the num ber of possible S ^ I transitions per time-step, 
and I{t) is the mean number of infected nodes. 


Empirically relevant networks are s par se an d t ransi- 
tion rates are small, so typically Q{t) E{t) ^ M{t). 
(The first inequality is a consequence of transition rates 
being small compared to 1/At; the second inequality fol¬ 
lows by noting that /(t) ^ A ^ ^(^)-) This means that 
the performance of the rejection sampling algorithm is 


limited by the rejection sampling step [second term of 
Eq. (13)], while the performance of the temporal Gille¬ 
spie algorithm is limited by the iteration over the set of 
contacts in order to update VL(t) [first term in Eq. (14)]; 
this explains why the difference in performance decreases 
with the mean instantaneous degree of the network. This 
also hints at how we may improve the speed of the tem¬ 
poral Gillespie algorithm: by rendering the identification 
of relevant contacts during each time-step faster. One 
such approach which may be applied to processes with 
an absorbing state (e.g. an IZ state) is explored below. 


Improving performance by removing obsolete contacts. 
Empirical networks describing human contact differ from 
simulated networks in a number of ways. Eor example, 
their structure and dynamics are more complex [25l[46]- 
l49] but perhaps most importantly in the perspective of 
optimizing simulations, they are of finite length. One is 
often interested in long-time behavior or slowly evolving 
processes compared to the length of available data. To 
overcome this limitation, one usually loops over the data 
set. This means that if a node enters an inactive ab¬ 
sorbing state such as the recovered (IZ) state in the SIR 
model, one may remove all following contacts to this node 
from the data, thus reducing the number of contacts that 
one must go through during the following loop. Eurther- 
more, since the X ^ IZ transition is independent of the 
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FIG. 5. Comparison of the speed of the temporal Gillespie and rejection sampling algorithms on empirical 
time-varying networks. Ratio ©rs/©tga between the time per realization of a single simulation using rejection sampling 
and the temporal Gillespie algorithm on empirical face-to-face contact networks in different social settings (Table |^: (A) for a 
SIR process, without (TGA) and with (TGA+GR) contact removal; (B) for a SIS process. Simulations were performed with 
/3 = 1 000/i for the workplace and /3 = lOO/x for the other networks. (Details on how simulations were performed are found in 
Methods.) 


network, one may also remove all contacts between two 
infected nodes. 

Pseudocode for an algorithm that removes obso¬ 
lete contacts is given in Methods and C++ code 
can be found at https://github.com/CLVestergaard/ 
TemporalGillespieAlgorithm, Figure [^A) compares 
the speed gain of the temporal Gillespie algorithm rela¬ 
tive to rejection sampling with and without contact re¬ 
moval for simulations of a constant-rate SIR process on 
empirical networks of face-to-face contacts (Table [^. De¬ 
pending on the parameters of the simulated process, re¬ 
moving obsolete contacts may induce both a significant 
gain or loss in speed; for processes that are fast compared 
to the length of the data set, the data is not repeated or 
only repeated few times during a simulation and the ad¬ 
ditional overhead involved in identifying and removing 
the obsolete contacts renders the algorithm slower; for 
slow processes the data is looped many times and re¬ 
moving the obsolete contacts makes the algorithm faster. 
Figure A) suggests an empirically determined rule-of- 
thumb: if the slowest time-scale of the simulated process 
(here ^ l//i) is longer than the length of the data, T, 
removing obsolete contacts pays off, if it is shorter, one 
should not remove obsolete contacts. 

Slow network dynamics. For time-varying networks 
of face-to-face contacts, which are relevant for simulat¬ 
ing epidemic spreading in a population, network dynam¬ 


ics are typically much faster than the time-scales of the 
dynamical process that is simulated (compare the 20 s 
time-resolution of the empirical data of Table [I] to typi¬ 
cal 1/^5 ^ 1 hour and l//i ^ 1 week for flu-like diseases). 
In the opposite case, i.e., if the network evolves much 
slower than the dynamical process, the temporal Gille¬ 
spie algorithm simply works like a static Gillespie algo¬ 
rithm in-between changes in the network structure while 
taking the changes changes into account exactly when 
they occur. The performance of the temporal Gillespie 
algorithm then approaches that of a st atic G ille spie al- 
gorithm in this case. Note that since Q{t) E{t) in 
this limit, the second term dominates in Eq. (14), which 
means that the speed of the algorithm is limited by the 
selection of waiting times and transitions that take place, 
and care should be taken to optimize these steps, e.g., by 
organizing the transition processes in a heap or a prior¬ 
ity queue m- Note finally that to obtain reliable re¬ 
sults using a rejection sampling algorithm one must use 
a time-step for simulations Atps which is much smaller 
than the time-step At of network change. Thus the ex¬ 
pected time complexity of rejection sampling scales with 
At/Atps ^simu ^ ^simu this casc. 
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F. Non-Markovian processes 


For real-world contagion processes, transition rates are 
typically not constant but in general depend on the his¬ 
tory of the process [301 [3l]. Such processes are termed 
non-Markovian. The survival probability for a single 
non-Markovian transition process taking place on a time- 
varying network is given by: 




= exp 





(15) 

Here is a filtration for the process m, i.e., all infor¬ 

mation relevant to the transition process available up to 
and including time t; typically, will be its starting 

time and relevant contacts that have taken place since. 
As in Sec. II D[ t* is the time of the last transition and 
t** = t* + r is the time of the next. [Note that since Xm 
now depends explicitly on t, we may absorb I^ in A^; 
however, to underscore the analogy with the Poissonian 
case, we keep the factor I^ explicitly in Eq. (15).] 


We use again that the transition processes are inde¬ 
pendent, to write the waiting time survival probability: 


S = exp 



with 




m^Vl{t) 


A (t; J^t) dt 




and where Tt is the union over for m G fl. 


(16) 


(17) 


For a static network, Eq. (16) reduces to the result 


found in na Eq. (7)]. This can be seen by noting that 
M(t) = M and = ft are then constant, and thus 

that A„(t; J-y)) = = 

d{ln[l/5„(t; and = S^it + 


yielding directly Eq. (7) of [iTl- 


-.(m)x 


As in the Poissonian case (Sec. IID) we define the nor¬ 
malized waiting time, r', as 


= K{t-Ft)dt . (18) 

This gives us the same simple forms as above for the 
survival function of the normalized waiting time, r'. 


S'(t') = exp(-T') , (19) 


and the probability that m is the transition that takes 
place at t = t**, 

Xm 

— ^m{t) \(^t’ J- ) * 

Until now our approach and results are entirely equiv¬ 
alent to the Poissonian case considered above. However, 


since A^(t) in general depend continuously on time, the 
transition time t** is not simply found by linear interpo¬ 
lation as in Eq. ( [T^ . Instead, one would need to solve the 
implicit equation L(t**;t*) = r' numerically to find t** 
exactly. To keep things simple and speed up calculations, 
we may approximate A(t) as constant over a time-step. 
This assumes that AA{t)At 1, where AA(t) is the 
change of A(t) during a single time-step. It is a more le¬ 
nient assumption than the assumption that A{t)At 1 
which rejection sampling relies on, as can be seen by not¬ 
ing that in general AA{t)/A{t) 1. The same assump¬ 
tion also lets us calculate L(tn+i; U) as in the Poissonian 
case: 


L(f„+i; r, 7)) = (f„.+i - f*)A(f*) +At ^(^d ’ 

i=n*+l 

( 21 ) 

and the time, t**, at which the next transition takes 
place: 


I** _ I 

i — ir) 


t' - 

A{tn**;Tt) 


( 22 ) 


Using the above equations, we can now construct a tem¬ 
poral Gillespie algorithm for non-Markovian processes. 

This algorithm updates all Xm{t) that depend on time 
at each time-step, where for the Poissonian case we only 
had to initialize new proces ses, i .e., contact-dependent 
processes (type b and c. Sec. H A| ). This means the algo¬ 
rithm is only roughly a factor two faster than rejection 
sampling [compare dotted lines (e = 0) in Fig. [^. To 
speed up the algorithm, we may employ a first-order cu- 
mulant expansion of S{r;J^t) around r = 0, as proposed 
in nmsH] for static non-Markovian Gillespie algorithms. 
It consists in approximating A^(t; by the constant 

A^(t*; for t* < t < t** and gives a considerable 

speed increase of the algorithm [full lines (e ^ oo) in 
Fig.§. However, the approximation is only valid when 
M{t) ^ I [43], which is not always the case for conta¬ 
gion processes. Notably, at the beginning and end of an 
SIR process, and near the epidemic threshold for an SIS 
process, M is small and the approximation breaks down; 
the approximate algorithm for example overestimates the 
peak number of infected nodes in a SIR process with re¬ 
covery rates that increase over time [compare full black 
line U ^ oo) to the quasi-exact full red line (e = 0) in 
Fig. I^]. An intermediate approach, which works when 
the number of transition processes is small, but is not too 
slow to be of practical relevance, is needed. We propose 
one such approach below [44] . 

Efficient non-Markovian temporal Gillespie algorithm. 
As discussed above, we neither want to update all tran¬ 
sition rates at each time-step as this makes the temporal 
Gillespie algorithm slow, nor do we want to only update 
them when a transition event takes place as this makes 
the algorithm inaccurate. 

An intermediate approach is found by looking at the 
relevant physical time-scales of the transition processes: 
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FIG. 6. Comparison of the speed of the temporal Gillespie and the rejection sampling algorithms: non- 
Markovian SIR process. Ratio ©rs/©tga between the time per realization of a single simulation of an SIR process 
with Weibull distributed recover y tim es using rejection sampling and the temporal Gillespie al^rithm on activity driven net¬ 
works of different average degree k{t) [for empirical contact networks, k(t) ~ 0.004-0.07 (Table W]: (A) for networks consisting 
of = 100 nodes and (B) of = 1 000 nodes. The parameter e controls the accuracy of the temporal Gillespie algorithm: 
for e = 0, where Am(t) is approximated as constant over a single time-step, it is most accurate; for e ^ oo, where Am is 
approximated as constant between two consecutive transition events, it is the least accurate. Node activities were given by 


di = TIZi with Zi ~ z^ 


for 


G [0.03,1); a node formed two contacts each time it was active. The recovery rate of an infected 
node was given by Eq. (23) with 7 = 1.5. The infection rate was [3 = 10^/io for networks with k{t) = 0.002, /3 = 10^/io for 
networks with k{t) = 0.02, and [3 — lO/xo for networks with k{t) = 0.2. (See Methods for details on how simulations were 
performed.) 




FIG. 7. Comparison of the outcome of non-Markovian SIR processes for different values of the parameter e. 

(A) Average number of infected nodes, (/), as function of time for a SIR process with Weibull distributed recovery times. (B) 
Distribution of the numbers of recovered nodes after the infection has died out (i.e. when 7 = 0). For e = 0 the temporal 
Gillespie algorithm is quasi-exact (see Supplementary Fig.|^for comparison with rejection sampling); for e ^ 00, corresponding 
to a first-order cumulant expansion of around t = t* (see main text), it is least accurate. As e is decreased, both (7)(t) 

and p{R) rapidly approach the quasi-exact result obtained for e = 0. Simulations were performed on an activity-driven network 
consisting of A^ = 100 nodes with activities ai = ^z/10, where Zi ~ nodes’ recovery times followed Eq. (23) with 7 = 1.5 

and po = 10“^ Hz; the length of a time-step was At = 1 s and the infection rate [3 = lOO/xo = 10“^ Hz. 
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the average waiting time before they take place, 

If the time elapsed since we last updated A^(t) is small 
compared to we do not make a large error by 

treating it as constant over the interval; however, if the 
elapsed time is comparable to or larger than the 

error may be considerable. Thus, instead of updating 
Xm at each time-step, we may update it only after a time 
t > has elapsed since it was last updated. Here e 

controls the precision of the algorithm. 

Below, we use this approach to simulate a 
non-Mar kovian SIR process, where the recovery 
times of infected nodes follow a Weibull distri¬ 
bution (see Methods for an algorithm written in 
pseudocode and https://github.com/CLVestergaard/ 
TemporalGillespieAlgorithm for implementation in 
C++). The recovery rate of an infected node is here 
given by 

, (23) 

where /io sets the scale, is time when the node was 
infected, and 7 is a shape parameter of the distribution. 
For 7 = 1, we recover the constant-rate Poissonian case 
with /i(t; = /io- For realistic modeling of infections, 

7 > 1 ; here is zero at t = and grows with 

time. In this case, we thus update the recovery rates 
whenever the time elapsed since a transition 
last took place exceeds = r(I + I/ 7 )//io. 

The parameter e lets us control the precision of the 
non-Markovian temporal Gillespie algorithm: the smaller 
e is, the more precise the algorithm is, on the other hand, 
the larger e is, the faster the algorithm is (Fig[^. At 
e = 0 , the temporal Gillespie algorithm is maximally ac¬ 
curate, but also slowest, corresponding to the quasi-exact 
approximation that stays constant over a single 

time-step. Letting e ^ 00 corresponds to the first or¬ 
der cumulant expansion of m, and is the fastest, but 
least accurate. Intermediate e gives intermediate accu¬ 
racy and speed, and permits one to obtain the desired 
accuracy without sacrificing performance. In the case of 
the SIR process with Weibull-distributed recovery times, 
e = O.I gives an error of no more than a few percent 
(Figs.[^-[^ and[^ —which is usually acceptable as the 
discrepancy between model and reality can be expected 
to be larger—with an almost optimal computation time 
(Figs. and|^. 


III. DISCUSSION 

We have presented a fast temporal Gillespie algorithm 
for simulating stochastic processes on time-varying net¬ 
works. The temporal Gillespie algorithm is up to mul¬ 
tiple orders of magnitude faster than current algorithms 
for simulating stochastic processes on time-varying net¬ 
works. For Poisson (constant-rate) processes, where it is 
stochastically exact, its application is particularly sim¬ 
ple. The algorithm is also applicable to non-Markovian 




FIG. 8. Accuracy and speed of the non-Markovian 
temporal Gillespie algorithm as function of e. (A)-(D) 
Different measures of the difference in outcome of simulations 
between algorithms with e > 0 and e = 0 (quasi-exact). (A) 
Difference, Azmax = ^max(e) — ^max(0), in the peak average 
fraction of infected nodes, Zmax = Vmax/A. (B) Difference 
Atmax between the times at which this peak takes place, nor¬ 
malized by /io- (C) Difference, Afoo, in the average frac¬ 
tion of nodes affected by the infection—the average attack 
rate. (D) Kullback-Leibler divergence, KL[p(roo)], between 
the distributions of attack rates (E) Time per simulation of 
the process. Simulations were performed on an activity-driven 
network with N = 100 nodes and activities = ZijXt) with 
Zi ~ for Zi G [0.03,1); nodes’ recovery times followed 

Eq. (23) with 7 = 1.5, /io = 10“^ Hz, and At = I s. 
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processes, where a control parameter lets one choose the 
desired accuracy and performance in terms of simulation 
speed. We have shown how to apply it to compartmen- 
tal models of contagion in human contact networks. The 
scope of the temporal Gillespie algorithm is more gen¬ 
eral than this, however, and it may be applied e.g. to 
diffusion-like processes or systems for which a network 
description is not appropriate. 
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METHODS 

The following four subsections contain supporting in¬ 
formation to the manuscript: the first subsection lists 
notation used in the article (Notation); the second de¬ 
tails how Monte Carlo simulations were performed (De¬ 
tails on how Monte Carlo simulations were performed) 
the third gives pseudocode for application of the tem¬ 
poral Gillespie algorithm to specific contagion processes 
on time-varying networks (Algorithms for simulating spe¬ 
cific contagion models). Finally, in the fourth subsection 
we give pseudocode for further optimization of the al¬ 
gorithm for empirical networks by removal of obsolete 
contacts (Removing obsolete contacts for an SIR process 
on empirical networks). 


Notation 

Tables m and mil list the notation used in the 
manuscript. Table |n| gives notation pertaining to the 
temporal Gillespie algorithm, and Table [TTT| lists notation 
pertaining to time-varying networks and compartmental 
contagion processes. 


Details on how Monte Carlo simulations were 
performed 

All simulations for comparing the speed of algorithms 
were performed sequentially on a HP EliteBook Folio 
9470m with a dual-core (4 threads) Intel Core i7-3687U 
CPU @ 2.10 GHz. The system had 8 GB 1600 MHz 
DDR3 SDRAM and a 256 GB mSATA Solid State Drive. 
Code was compiled with TDM GCC 64 bit using g++ with 


the optimization option -02. Speedtests were also per¬ 
formed using -03 and -Ofast, but -02 gave the fastest 
results, both for rejection sampling and temporal Gille¬ 
spie algorithms. 

For SIR processes simulations were run until / = 0; for 
SIS processes simulations were run for 20/(//At) time- 
steps (as in Fig. or until / = 0, whichever happened 
first. Between 100 and 10 000 independent realizations 
were performed for each data point depending on fiAt 
(100 for low /iAt and 10 000 for high jiAt). For simula¬ 
tions on empirical contact data, data sets were looped if 
necessary. 


Algorithms for simulating specific contagion models 

We here give pseudocode for the application of the tem¬ 
poral Gillespie algorithm to three specific models: the 
first subsection treats the Poissonian SIR process, the 
second treats the Poissonian SIS process, and the third 
treats a non-Markovian SIR process with recovery times 
following a general distribution. 

We assume in the following that the time-varying net¬ 
work is represented by a list of lists of individual con¬ 
tacts taking place during each time-step. An individ¬ 
ual contact, termed contact, is represented by a tu¬ 
ple of nodes, i and j. The list contactLists [t] gives 
the contacts taking place during a single time-step, t, 
fort=0,lj- • . ,T_simulation-1, where T_simulation is 
the desired number of time-steps to simulate. The state 
of each node is given by the vector x, where the entry 
X [i] G {S, I j h} gives the state of node i. 

As one may always normalize time by the duration of a 
time-step. At, we have in the following set At = 1. Note 
that beta and mu in the code then corresponds to 
and /iAt, respectively. 

SIR process. The classical SIR model with constant 
infection and recovery rates is the simplest epidemic 
model where individuals can gain immunity. As dis¬ 
cussed in the main text, nodes may be in one of three 
states: susceptible (5), infectious (X), or recovered {IZ). 
Two different transition types let the nodes switch state: 
a spontaneous I ^ IZ transition which takes place with 
rate //, and a contact-dependent S ^ X transition which 
takes place with rate P upon contact with an infectious 
node. Pseudocode shows how the temporal Gillespie 
algorithm may be implemented for an SIR process on a 
time-varying contact network. 

SIS process. In the SIS model, nodes can be in one of 
two states: susceptible (S) or infectious (X). As for the 
SIR model, two different transition types let the nodes 
switch state: a spontaneous X ^ S transition which takes 
place with rate /i, and a contact-dependent S ^ X tran¬ 
sition which takes place with rate P upon contact with 
an infectious node. 

The SIS model is implemented in a manner very sim¬ 
ilar to the SIR model; an implementation can be found 
by using the code of Pseudocode with lines 07 and 
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TABLE IL Notation pertaining to the temporal Gillespie algorithm. 


Symbol 

Description 

First appearance(s) 

t 

Real time. 

Sec. 

IIA 


At 

Duration of a time-step. 

Sec. 

IIB 


n 

Time-step number. 

Sec. 

IID 


t"n 

Time at beginning of time-step n: tn = nAt. 

Sec. 

IID 


m 

Possible transition / transition process. 

Sec. 

IIA 


Am 

Transition rate for m. 

Sec. 

IIA 


Im(t) 

Eunction indicating if the transition m may take place at time t. 

Sec. 

IID 


n{t) 

Set of transition processes at time t. 

Secs. 

llA 

lIID 


M{t) 

Number of transition processes at time t. 

Secs. 

IIA 

lIID 


Q 

Set of total possible transitions between two consecutive transition events. 

Secs. 

llA 

lIID 


M 

Number of total possible transitions between two consecutive transition events. 

Secs. 

IIA 

lIID 


A, A(i) 

Cumulative transition rate (at time t): A(t) = X)mGO(t) 

Secs. 

IK 

UID 


L(t; t*) 

Integrated cumulative transition rate (from t* to t). 

Sec. 

IID 


T 

Waiting time between two consecutive transitions. 

Sec. 

lie 


S{t) 

Waiting time survival function. 

Sec. 

lie 



Times when the last/next transition took/takes place, respectively. 

Sec. 

IID 


n*, n** 

Time-steps during which the last/next transition took/takes place, respectively. 

Sec. 

IID 


T 

Normalized waiting time between two consecutive transition events. 

Sec. 

IID 


S{r') 

Normalized waiting time survival function. 

Sec. 

IID 


t' ~ Exp(I) 

T is exponentially distributed with unit rate. 

Sec. 

IID 

(Sec. lie I 

©RS 

Time per simulation for the rejection sampling algorithm. 

Sec. 

HE 


©TGA 

Time per simulation for the temporal Gillespie algorithm. 

Sec. 

HE 


0{x) 

Term of order x, i.e., O (x) — ax for a given constant a. 

Sec. 

HE 


-r{m) 

t 

Eiltration for the transition process m. 

Sec. 

HE 


Tt 

Union of all 

Sec. 

HE 



40 removed and line 37 replaced by x[m] = S. C++ 
code is found at https://github.com/CLVestergaard/ 
TemporalGillespieAlgorithm for both homogeneous 
and heterogeneous populations. 

Non-M arkov ian SIR process. We consider in the main 
text (Sec. ES) an SIR model with non-constant recov¬ 
ery rates; instead of being exponentially distributed (as 
in the constant-rate SIR model), the individual recovery 
times, are here Weibull distributed. 


(m) 


’ 7/^0 




7-1 


O-MOT 


(m) 


(24) 


As for the classical SIR model, nodes may be in one of 
three states: susceptible (5), infectious (X), or recovered 
{IZ). Two different transition types let the nodes switch 
state: a contact-dependent S ^ X transition with con¬ 
stant rate (3 upon contact with an infectious node, and 
a spontaneous X ^ IZ transition which takes place with 
rate given by Eq. (23). 


The implementation of the SIR model with non- 
exponentially distributed waiting times requires some ex¬ 
tension of the code for the constant-rate SIR model to ac¬ 
count for the heterogeneous and time-dependent recovery 
rates. To this end, we introduce the following variables: 


t_inf lists the times at which each node was infected (if 
applicable); t* is the exact time at which the last tran¬ 
sition took place; mu is a function of time that is called 
as mu(t-t_inf [m]) to return the instantaneous recov¬ 
ery rate of m at time t; mu_avg is the expected recovery 
time of an infected node and is used together with the 
precision control parameter epsilon in the approximate 
simulation scheme discussed in Sec. in FI Pseudocode [2] 
shows pseudocode for an implementation of such a SIR 
model with non-constant recovery rates. 


Removing obsolete contacts for an SIR process on 
empirical networks 


When simulations are carried out on data which are 
looped due to their finite length, the speed of the tem¬ 
poral Gillespie algorithm may be further increased for 
processes with an absorbing state, such as the SIR pro¬ 
cess by removing obsolete contacts to nodes that have 
entered such a state. Pseudocode shows pseudocode 
for removing obsolete contacts; its replaces lines II to 19 
of Pseudocode [T] 
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TABLE III. Notation pertaining to compart mental contagion models and time-varying networks. 


Symbol 

Description 

First appearance 

L 3 

Node. 

Sec. 

HA 


N 

Number of nodes in network. 

Sec. 

HA 



Contact taking place at time t between nodes i and j. 

Sec. 

HA 


kxit) 

Number of infected nodes in contact with i at time t 

Sec. 

HA 


k{t) 

Average degree (number of contacts per node) of network at time t. 

Fig. 

4 

Xi{t) 

Random variable specifying the state (compartment) of node i at time t. 

Sec. 

HA 


xe{Xi,X2...Xq} 

Possible node states (compartments). 

Sec. 

HA 


Xj, 

Number of nodes in state Ap. 

Sec. 

HA 


<S, X, 

Possible node states in SIS and SIR models of epidemic spreading. 

Sec. 

HA 


S', /, R 

Number of nodes in each of the states <S, X, and 77., respectively. 

Sec. 

HA 


p 

Rate of <S ^ X transition of a susceptible node in contact with an infectious node. 

Sec. 

HA 



Rate of spontaneous X ^ IZ or X ^ S transition of an infectious node. 

Sec. 

HA 


m 

Mean number of contacts during a single time-step. 

Sec. 

HE 


M{t) 

Mean number of transition processes per single time-step. 

Sec. 

HE 


Q(t) 

Mean number of transitions that take place per time-step. 

Sec. 

HE 



Mean number of S-X contacts during a single time-step. 

Sec. 

HE 


m 

Mean number of infectious nodes. 

Sec. 

HE 


T 

Length of a data set describing a time-varying network (in time). 

Sec. 

HE 


"^simu 

Number of time-steps that are simulated during a single realization. 

Sec. 

HE 


AtRs 

Time-step used for rejection sampling when Am At are large, Atns < At. 

Sec. 

HE 


Mo 

Scale parameter of the Weibull distribution. 

Sec. 

HE 


7 

Shape parameter of the Weibull distribution. 

Sec. 

HE 


tim) 

Starting time for transition process m (e.g. the time when a node was infected). 

Sec. 

HE 
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//Initialize: 

01 FOR i=l,...,N 

02 x[i] = S //set node states to S 
03 ENDFOR 

04 x[root] = I //set state of root node to I 
05 m_I = [root] //list of infected nodes 
06 N_I = 1 //number of infected nodes 
07 N_R = 0 //number of recovered nodes 
08 Mu = mu //cumulative recovery rate 
09 tau = randexp(l) //draw tau ~ Exp(l) 

//Run through the time-steps: 

10 FOR t=0,l,...,T_simulation-l 

//Update list of possible S->I transitions: 

11 CLEAR m_SI //S nodes in contact with I nodes 

12 FOR contact in contactLists[t] 

13 (i,j) = contact 

14 IF (x[i] ,x[j] )==(S,I) 

15 APPEND i to m_SI 

16 ELSE IF (x[i] ,x[j])==(I,S) 

17 APPEND j to m_SI 

18 ENDIF 

19 ENDFOR 

20 M_SI = length of m_si 

21 Beta = beta*M_SI //cumulative infection rate 

22 Lambda = Mu+Beta //cumulative transition rate 


//Check if a transition takes place: 

23 IF Lambda<tau //no transition 

24 tau -= Lambda 

25 ELSE //at least one transition 

26 xi = 1. //remaining fraction of time-step 

27 WHILE xi*Lambda>=tau 

28 DRAW z uniformly from [0,Lambda) 

29 IF z<Beta //S->I transition 

30 DRAW m at random from m_SI 

31 x[m] = I 

32 APPEND m to m_I 

33 N_I += 1 

34 Mu += mu 

35 ELSE //I->R transition 

36 DRAW m at random from m_I 

37 x[m] = R 

38 REMOVE m from m_I 

39 N_I -= 1 

40 N_R += 1 

41 Mu -= mu 

42 ENDIF 

43 xi -= tau/Lambda //update remaining fraction 

//Update list of S->I transitions and rates: 

44 REDO lines 11-22 

45 tau = randexp(l) //draw new tau 

46 ENDWHILE 

47 ENDIF 

//Read out the desired quantities: 

48 WRITE N_I, N_R, ... 

49 ENDFOR 


Pseudocode 1. Pseudocode for an SIR process with constant and homogeneous transition rates. C++ code for 
homogeneous and heterogeneous populations is found at https://github.com/CLVestergaard/TemporalGillespieAlgorithm 
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//Initialize: 

01 FOR i=l,...,N 

02 x[i] = S //set nodes states to S 
03 ENDFOR 

04 x[root] = I //set state of root node to I 
05 t_inf[root] = 0 //time of infection = 0 
06 m_I = [root] //list of infected nodes 
07 mus = [mu(0)] //list of recovery rates 
08 Mu = mu(0) //cumulative recovery rate 
09 N_I = 1 //number of infected nodes 

10 N_R = 0 //number of recovered nodes 

11 tau = randexp(l) //draw tau ~ Exp(l) 

//Run through the time-steps: 

12 FOR t=0,l,...,T_simulation-l 

//Update mus if t-t*>=epsilon*mu_avg: 

13 IF t-t*>=epsilon*mu_avg 

14 CLEAR mus 

15 FOR m in m_I 

16 APPEND mu(t-t_inf[m]) to mus 

17 ENDFOR 

18 Mu = sum of mus 

19 ENDIF 

//Update list of possible S->I transitions: 

20 CLEAR m_SI //S nodes in contact with I nodes 

21 FOR contact in contactLists [t] 

22 (i,j) = contact 

23 IF (x[i] ,x[j])==(S,I) 

24 APPEND i to m_SI 

25 ELSE IF (x[i] ,x[j] )==(!, S) 

26 APPEND j to m_SI 

27 ENDIF 

28 ENDFOR 

29 M_SI = length of m_si 

30 Beta = beta*M_SI //cumulative infection rate 

31 Lambda = Mu+Beta //cumulative transition rate 


//Check if transition takes place: 

32 IF Lambda<tau //no transition 

33 tau -= Lambda 

34 ELSE //at least one transition 

35 xi = 1. //remaining fraction of time-step 

36 t* = t //for calculating transition times 

37 WHILE xi*Lambda>=tau 

38 t* += tau/Lambda //transition time 

39 DRAW z uniformly from [0,Lambda) 

40 IF z<Beta //S->I transition 

41 DRAW m at random from m_SI 

42 X [m] = I 

43 t_inf[m] = t* 

44 APPEND m to m_I 

45 N_I += 1 

46 ELSE //I->R transition 

47 DRAW m from m_I with weight mus[m] 

48 x[m] = R 

49 REMOVE m from m_I 

50 N_I -= 1 

51 N_R += 1 

52 ENDIF 

53 xi -= tau/Lambda //update remaining fraction 

//Update mus: 

54 CLEAR mus 

55 FOR m in m_I 

56 APPEND mu(t*-t_inf[m]) to mus 

57 ENDFOR 

58 Mu = sum of mus 

//Update list of S->I transitions and rates: 

59 REDO lines 20-31 

60 tau = randexp(l) //draw new tau 

61 ENDWHILE 

62 ENDIF 

//Read out the desired quantities: 

63 WRITE N_I, N_R, ... 

64 ENDFOR 


Pseudocode 2. Pseudocode for a non-Markovian SIR process with non-constant recovery rates. The function mu 
returns the instantaneous recovery rate as function of (t — t*); for Weibull distributed recovery times, mu is given by Eq. (23). 
C++ code is found at https://github.com/CLVestergaard/TemporalGillespieAlgorithm 
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01 CLEAR m_SI //S nodes in contact with I nodes 
02 FOR contact in contactLists[t] 

03 (i,j) = contact 

04 IF x[i]==S 

05 IF x[j]==I 

06 APPEND i to m_SI 

07 ELSE IF x[j]==R //remove if x[j]==R 

08 REMOVE contact from contactLists[t] 

09 ENDIF 

10 ELSE IF x[i]==I 

11 IF x[j]==S 

12 APPEND j to m_SI 

13 ELSE //remove if (x[i],x[j])==I or x[i]==R 

14 REMOVE contact from contactLists[t] 

15 ENDIF 

16 ELSE //remove if x[i]==R 

17 REMOVE contact from contactLists[t] 

18 ENDIF 

19 ENDFOR 

Pseudocode 3. Pseudocode for counting possible S ^ T transitions with removal of outdated contacts. C++ code 
is found at https://github.com/CLVestergaard/TemporaIGiIlespieAIgorithm 


Supplementary Figures 
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Supplementary FIG. 1. Numerical results from temporal Gillespie and rejection sampling algorithms for 
contagion processes taking place on empirical networks. (A)-(D) for a SIR process and (E)-(H) a SIS process. 
(A),(B),(E), and (F) for /3At = 10“^ and /xAt = 10“^; (C),(D),(G), and (H) for /3At = 10“^ and /xAt = 10“^. (A),(C) Mean 
number of nodes in each state of the SIR model as function of time. (B),(D) Distribution of hnal epidemic size (number of 
recovered nodes when / = 0) in the SIR model. (E),(G) Mean number of nodes in each state of the SIS model as function of 
time. (F),(H) Distribution of the number of infected nodes in the stationary state (t ^ oo) of the SIS model. All simulations 
were performed 1 000 000 times with the root node chosen at random on a face-to-face contact network recorded in a high school 
(Table l^l. 
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Supplementary FIG. 2. Comparison of numerical results from temporal Gillespie and rejection sampling algo¬ 
rithms for high transition probability per time-step. (A)-(D) for a SIR process and (E)-(H) a SIS process. (A),(B),(E), 
and (E) for /3At = 10“^ and /xAt = 10“^; (C),(D),(G), and (H) for /3At = 1 and /xAt = 10“^. (A),(C) Mean number of nodes 
in each state of the SIR model as function of time. (B),(D) Distribution of final epidemic size (number of recovered nodes 
when / = 0) in the SIR model. (E),(G) Mean number of nodes in each state of the SIS model as function of time. (E),(H) 
Distribution of the number of infected nodes in the stationary state (t oo) of the SIS model. All simulations were performed 
1 000 000 times with the root node chosen at random on an activity driven network consisting of A = 100 nodes, with activities 
di = where 77 = 0.1 and Zi ~ for zi G [0.03,1), and a node formed two contacts when active. 
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Supplementary FIG. 3. Comparison of numerical results from temporal Gillespie and rejection sampling al¬ 
gorithms for a non-Markovian SIR process. (a),(c) Mean number of nodes in each state as function of time in the SIR 
model with Weibull distributed recovery times (Sec. VII A); the parameter controlling the precision of the temporal Gillespie 
algorithm was set to e = 0 (quasi-exact). (b),(d) Distribution of final epidemic size (number of recovered nodes when / = 0). 
(a),(b) /3At = 10“^ and jiAt — 10“^; (c),(d) /3At = 10“^ and jiAt — 10“^. The outcome of the rejection sampling algorithm 
approaches that of the temporal Gillespie algorithm as /3At and /xAt become smaller. All simulations were performed 100 000 
times with the root node chosen at random on an activity driven network consisting oi N — 100 nodes, with activities ai — r]Zi, 
where 77 = 0.1 and Zi ~ for Zi G [0.03,1), and a node formed two contacts when active. Nodes’ recovery times followed 

Eq. (20) with 7 = 1.5 and the length of a time-step was At = 1 s. 















































